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Wang-Landau Simulation for the Quasi-One-Dimensional Ising Model 
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We revisit the nature of the quasi-one-dimensional Ising model on the basis of Wang-Landau 
simulation. We introduce the density of states in the two-dimensional energy space correspond- 
ing to the intra- and interchain directions. We then analyze the interchain coupling dependence 
of specific heat of the anistropic two-dimensional Ising model in the context of the density of 
states, and further discuss the size dependence of peak temperature. We also discuss the feature 
of the phase transition in the three-dimensional case. 
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1. Introduction 

The effect of the interchain coupUng in the quasi-one- 
dimensional(QlD) spin system has recently attracted 
much attention. In general, the ID spin system shows 
no phase transition at a finite-temperature. In the QID 
system, however, the weak interchain coupling induces 
a finite-temperature phase transition. In fact, the 3D 
long-range ordering for the QID system has been a 
long-standing issue. ^ The recent experimental develop- 
ments enable the precise investigation of 3D ordering in 
a wide variety of QID systems, such as coupled S = 1/2 
Heisenberg chains.^ Moreover, very recently, it has been 
shown that the QID Ising-like XXZ antiferromagnet 
BaCu2V208 exhibits an exotic incommensurate spin or- 
der in a magnetic field, ^■'^ in which the Ising anisotropy 
plays an essential role. 

Motivated by the experimental results above, we reex- 
amine the phase transition of the QID Ising system: 
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where S S ±1 is the Ising spin variable, < i,j >„ indi- 
cates spin pairs along the chain direction, and < i,j >±_ 
means the pairs perpendicular to the chain. Thus, J and 
J' respectively represent the coupling constants for the 
intra- and inter-chain interactions. Of course the univer- 
sality of the Ising transition due to the Z2 symmetry 
breaking itself is independent of the interchain interac- 
tion. However, quantitative details of the interchain de- 
pendence of the phase transition still involve an interest- 
ing problem, which is essential to resolve experimental 
results. Recently, the universal reduction in the effective 
coordination number has also been reported for the QID 
Ising model. ^ The aim of this paper is to understand the 
role of interchain coupling in the context of the energy 
density of states (DOS) based on Wang-Landau simula- 
tion.""^^' ^'^ The size dependence of the transition temper- 
ature of the QID Ising system is also discussed. 

For the quantitative analysis of the phase transition 
of the QID system, recall that the energy scale of inter- 
chain coupling is much smaller than that of intrachain 
coupling, and thus the transition temperature becomes 
very low; the conventional metropolis Monte Carlo simu- 



lation based on a local spin flip often fails in relaxation to 
the appropriate equilibrium state. Recently, an cfhcient 
cluster algorithm has been proposed for the QID sys- 
tem. -'^'^ In this paper, however, we employ such a gener- 
alized ensemble method as multicanonical simulation. 
In particular, the Wang-Landau simulation^^' enables 
us to estimate DOS efficiently through a random walk in 
the energy space, and to avoid trapping in a metastable 
state. Then we can resolve the contribution of typical 
configurations at low temperatures, which provides an 
essential viewpoint of the low-temperature behavior of 
the QID system. 

This paper is organized as follows. In the next section, 
we briefly explain details of the Wang-Landau simulation 
for the QID system. In particular, we introduce the DOS 
of two-dimensional energy space for the intra- and inter- 
chain directions. In §3, we discuss the phase transition 
in the 2D case in the context of DOS and then analyze 
the interchain interaction dependence of the phase tran- 
sition. For 3D case, we also discuss the nature of the 
phase transition. In §4, the summary and discussions are 
given. 

2. Simulation Details 

The Wang-Landau simulation is based on a random 
walk in the energy space without trapping metastable 
state and enables us to estimate DOS. For the QID 
system, however, the energy scale of interchain cou- 
pling is fairly different from that of intrachain coupling. 
For the purpose of treating such a highly anisotropic 
system more efficiently, we further introduce the two- 
dimensional energy space defined by 
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where e,, and e± respectively denotes the (dimensionless) 
unbiased energies for the intra- and inter-chain direc- 
tions. Then the total energy is given hy E = —Je„ — J'e±. 
The Wang-Landau simulation itself is performed for this 
two-dimensional space of the spatially isotropic Ising 
model and then DOS g{e„,e±) is obtained in (eii,ej^) 
space. The expectation value for various J' can be ob- 
tained through reweighting. Here, it should be noted that 
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the multicanonical simulation for the two-dimensional 
parameter space was successfully applied to such com- 
plex systems as the polymer system^"^ and the spin glass 
system. -'^^ Moreover, the Wang-Landau sampling drasti- 
cally enhances the accessibility to the multi-dimensional 
parameter space in various situations. An interesting 
point in the present case may be that the dimensionality 
of the parameter space is directly related to the spatial 
dimension. 

The detailed conditions for the Wang-Landau simu- 
lation are given as the follows. At the start of simula- 
tion, DOS is unknown, so it is simply set at g{e,„ e±) = 
1 for all possible (eii,ej^) . Then we begin a random 
walk in {e„,e±) space with a probability proportional to 
l/g{e,„e±) . The transition probability from (e||,ej^) to 
(e„ , e_^) is 



Prob((e||, e_i_) (e,, , e_|_)) = min 



g(eipe_L) ^ 
fi'(eM,el)' 



, (3) 



and g{e„,e±) is iteratively updated by the modification 
factor / as 



ln.g(e„,ej_) ^ ln.g(e„,e^) +ln/, 



(4) 



when the state is visited. At the same time, the his- 
togram is updated like H{e„,e±) — > H{e„, e±) + 1. When 
the histogram becomes "flat" , the modification factor is 
reduced, In/ ^ (ln/)/2, and we reset the histogram 
to zero. Then we perform a random walk again. To 
check the flatness of the histogram, wc use the criterion 

(-ffmax - -f^min)/ (-f^max + -ffmin) < 0.1 ~ 0.3, whcrC ffmax 

and Hjnin are the maximum and minimum histogram 

counts, respectively."'^^ We end the simulation when the 
modification factor is smaller than a predefined value (we 
set In f final = 10~^). The initial value of the modification 
factor is In/o = 1. The update of g{e„, ej_) and iJ(eii, e^) 
is performed every N spin flips, where N is the number 
of spins in the system.^'' 

3. Results 

3.1 2D Ising model 

Let us first consider the 2D Ising model on a L x L 
square lattice, the exact solution of which is well known^^ 
and very useful for verifying the simulation results. The 



Hamiltonian of the 2D Ising model is written as 

L L 

- J ^ SijSi+ij - J' ^ SijSij+i, (5) 

where i and j are the indexes of the intra- and inter- 
chain directions, respectively. The unbiased energies for 
the intra- and inter-chain directions are explicitly given 
by e,i = Y^tj St,jS,+i,j and e± = Yld,] Si,jSi^j+i, re- 
spectively. Since g{e»,eA_) = g{-e»,e±) = g{e,„-e±) = 
g{—e„, —e±) holds, it is sufficient to perform simulation 
in the region e„ ^ 0, ^ 0. The system sizes are 
L = 10, 20, 30, 40, and 50. The maximum histogram 
count per stage is .ffmax = 3024, and the maximum 
Monte Carlo steps per stage is 1.7 x 10^ for the L = 10 
system. Then the total number of stages is 21, and the to- 
tal CPU time is 2 minutes with a 2.66GHz Core2Duo pro- 
cessor. For L = 30, i?max = 7981, the maximum Monte 
Carlo steps per stage is 3.3 x 10^^. The total number of 
stages is also 21 and the total CPU time is 55 hours. In 
Fig. 1, we show the typical result of DOS g{e,„e±) for 
L = 10. 

On the basis of DOS g{e„,e±), we calculate specific 
heat; Figure 2 shows the size dependence of the specific 
heat for J'/ J = 0.025. According to the exact solution 
of the 2D Ising model, the transition temperature for 
J' /J 0.025 is given as Tc/J = 0.6221 ■ • • . The resuh 
clearly shows that the peak of C corresponding to the 
critical divergence gradually develops for L = 50 in the 
vicinity of Tc. In addition to the critical point, wc can also 
see a small peak in the low-temperature region T/J ~ 
0.2. As L increases, the peak temperature of this small 
peak shifts to a higher temperature side and the peak 
height itself decreases rapidly. 

In order to see the origin of the low-temperature peak, 
we show the DOS g{e„,e±) in the low-energy region in 
Fig. 3. Note that the scale of e,, is much smaller than that 
of e±; Since J' <^ J for the QID system, the range of the 
horizontal axis in Fig. 3 is adjusted by the ratio: en/e^ ~ 
J'/ J = 1/40. The ground state energy is E'g = -{J + 
J')L X L and its configuration is illustrated as Fig. 3(a), 
which is located at (e,,, e^) = (100, 100). As a low-energy 
excitation, the single-spin flipped state given by Fig. 3(b) 
is usually considered, whose energy is given as £^((,) = 
Eg + 4(J -|- J'). For the QID system, however, another 
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Fig. 3. Two-dimensional DOS of the 2D Ising model of L = 10 
near the groundstate(88 ^ e,, ^ 100,0 ^ ex ^ 100). The figures 
in the right panel indicate the configurations for (a) the ground 
state, (b) single-spin flipped state, and (c) chain flipped state. 




0.23 



Fig. 4. Energy distribution function for J'/ J = 0.025 with 
L = 10. The solid squares indicate the distribution function for 
T/J = 0.23, which corresponds to the low temperature peak 
of the specific heat. The open circles indicate that for the high 
temperature peak (T/J = 0.91). Each distribution function is 
normalized so that the maximum value corresponds to unity. 



important excitation we should discuss is "chain flipped 
excitation" , a typical example of which is depicted in 
Fig. 3(c), and its energy is given as = Eg + AJ'L. 
The DOS of the chain flipped configuration is located at 
e_L = 60, 20 • • • on the edge of e,, = 0. In Fig. 1, we can 
also confirm that the DOS of these configurations at the 
edges deviate from the "bulk" value in the (e,,, e±) plane. 
Moreover, note that the "gap" in DOS at every e± = 20 
in Fig. 3 also originates from the chain structure of the 
lattice. Since L < J/J', we can see that the chain flipped 
excitation has a lower energy than the single spin flipped 
state, but it becomes the predominant excitation at a low 
temperature. As L increases beyond J/ J', the energy of 
the chain flipped excitations shifts to the higher-energy 
region, so that the contribution of such configurations 
decreases gradually. Thus the low-temperature peak of 
specific heat in Fig. 2 can be well described by the chain 
flipped excitations, which is peculiar to the QID system. 

In order to see the weight of each energy state in the 
equilibrium, we calculate the energy distribution func- 
tion 

P (e,„ e^) = 5(e,„ e±) exp [(Je„ + J'e^)/T] (6) 



for J'/ J = 0.025, which is shown in Fig. 4. Note that 
T/J = 0.23 is the temperature of the low-temperature 
peak of specific heat and T/ J = 0.91 corresponds to the 
high-temperature peak of L = 10. In the figure, the pre- 
dominant contribution at T/J ~ 0.23 clearly comes from 
the states at the edge of e± = 100, which implies that the 
chain fiipped state is essential for the small peak; For a 
relatively small system size [L < J/ J'), the energy of the 
excitations actually satisfies 4( J + J') > AJ'L. On the 
other hand, P{e„, e±) for T/ J = 0.91 shows a Gaussian- 
like shape at approximately (eii,e_L) (90,20), where 
the chain fiipped state gives only a minor contribution 
in DOS. 

As mentioned above, the predominant contribution 
to the low-temperature peak(T/J = 0.23) is the chain 
fiipped states near the ground state. This implies that 
the polarization of the spins in the same chain is basi- 
cally frozen and the aligned spins in the chain can be- 
have as a single spin, which forms an effective ID spin 
chain through the weak interchain coupling LJ' in the 
interchain direction. Thus, we can see that the fluctu- 
ation in the interchain direction is predominant for the 
low-temperature peak, while for the high temperature 
peak, the fluctuations in both the intra- and inter-chain 
directions give the signiflcant contributions. Of course, 
the low temperature peak is basically a flnite-size effect 
and vanishes in the bulk limit. However, the region where 
the finite-size effect can clearly appear is up to L ~ J/J', 
which is a certain large number for the QID system. This 
implies that the true critical divergence of specific heat 
is eventually masked by the analytic contribution of the 
low-temperature peak, up to L ~ J / J' . Thus, the finite- 
size scaling analysis based on the data L < J / J' should 
be performed carefully. 

In order to extract the proper critical behavior for 
L < J / J' , we examine the decomposition of specific 
heat into three parts: the fluctuation along the chain, 
C,,; the fluctuation in the interchain direction, C^_; and 
cross term of the intra- and inter-chain directions, Ccross- 
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i{E^)- {E)'^)/NT^ 
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The intra- and inter-chain fluctuations basically reflect 
the ID nature of the spin fluctuations. Thus, for C, C,,, 
and Cj_, we have to see the critical fluctuation in the 
background of the predominant ID behavior. On the 
other hand, the cross term Ccross exhibits no such ID 
behavior and thus we can expect more direct observa- 
tion of the critical fluctuation. 

In Fig. 5, we show C±_ and Cdoss for J' / J = 0.025(C|| is 
not presented here). In the flgure, Cj_ shows a Schottky- 
like peak for small system sizes [L < 20). As L increases, 
the peak position shifts to the high-temperature side and 
the peak height itself rapidly decreases. This behavior is 
consistent with the fact that the interchain fluctuation is 
the predominant contribution for L < J / J' . Indeed, we 
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Fig. 5. Decomposed specific heats for J' j J = 0.025: (a) the en- 
ergy fiuctuation in the interchain direction, C^_, and the cross 
term of the chain and interchain directions, Ccross- 
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where A is a nominivcrsal constant . First wo discuss the 
results for J'/ J = 0.1, which are illustrated in Fig. 6(a). 
In the figure, we can see that C, C,,, and C_l gradually 
approach Tc, for which no scaling behavior is observed. 
However, the cross term of the specific heat Ccross well 
satisfies eq. (11) within a relatively small system size 
(1/L < 0.06), suggesting that Ccross is more effective for 
capturing the critical behavior than C. Another interest- 
ing feature of the QID system is that Tc{L) approaches 
Tc from the under side of T^, namely, A > 0, for suffi- 
ciently large L. This behavior is contrasted to ^ < in 
the isotropic case where the peak temperatures of all the 
Cs monotonically approach Tc from T > Tc- 

The peak temperatures for J'/ J = 0.025, which are 
shown in Fig. 6(b), demonstrate a more typical size de- 
pendence of the QID system; The peak of C,, monotoni- 
cally decreases from the upper side of Tc, while C± clearly 
exhibits the crossover behavior. For small L{< 0.04), the 
peak position of Cj_ originates from the chain flip config- 
uration. However, we can see that it rapidly crossovers 
to that of the critical behavior at 1/L ~ 0.04. On the 
other hand, Ccross seems to approach Tc smoothly, sug- 
gesting that Ccross is more suitable for finite size analysis 
of the critical behavior. For J'/ J = 0.025, however, note 
that system size may be still insufficient for the precise 
verification of the critical exponent v. 

3.2 3D Ising model 

Let us discuss the 3D Ising model along the same line 
of argument as that of the 2D case. The Hamiltonian is 
written as 

L 

'H = —J Sij^kSij,k+i 



Fig. 6. Size dependences of the peaJc temperatures for C, Cii, Cx, 
and Ccross: (a) J' /J = 0.1 and (b) J' /J = 0.025. The soUd circles 
at the vertical axis indicate the exact transition temperatures. 



have verified that the low temperature peak of i = 10 

can be well fitted by the specific heat of the ID Ising 
chain of the effective coupling LJ' with L = 10. In ad- 
dition to the low-temperature peak, a broad peak also 
emerges at approximately T/J ^ 0.6, as L increases; 
this peak corresponds to the critical divergence in the 
bulk limit. Thus, the crossover of C± from the effective 
ID Ising model behavior to the 2D Ising model clearly 
appears at L ~ J/J'. On the other hand, the cross term 
Ccross shows the divergence behavior only near the cor- 
rect critical temperature Tc = 0.622 • • • . This suggests 
that Ccross captures the critical behavior more effectively 
than the total specific heat C. 

Wc further analyze the size dependence of the peak 
temperatures of the decomposed specific heats C, d, 
C_L, and Ccross- Let us write the peak temperature for 
the system size L as Tc{L). Then, in the critical regime, 
peak temperature is expected to follow the size scaling 



Tc{L) -Tc = AL- 



(11) 



— J' '^^[Si,j,kSi-\.ij,k + Si^j^kSi,j+i,k\^ (12) 

where k is assumed to run in the chain direction. In ac- 
tual computations, the maximum histogram count per 
stage is ifmax = 6024, and the maximum Monte Carlo 
steps per stage is 6.5 x 10^ for L = 6 system. The to- 
tal number of stages is 21 and the total CPU time is 50 
minutes. For L = 10, J/max = 6748, the maximum Monte 
Carlo steps per stage is 7.5 x 10^^, and the total CPU 
time is about 6 days. Here, we note that, for the 3D Ising 
model, Wang-Landau simulation in the 2D energy space 
(2) sometimes does not achieve good convergence near 
the edges of the 2D energy space. In such a case, we have 
additionally performed Wang-Landau simulation for the 
conventional ID energy space to obtain specific heats. 

Figure 7(a) shows the size dependence of the specific 
heat C for J'/J = 0.025 up to T = 18. In the figure, 
we can see the broad maximum of C of i = 6 at ap- 
proximately T/J ^ 1.0. At the same time, a small peak 
emerges in the low-temperature region, refiecting the ID 
nature of the system. As L increases, this small peak 
rapidly merges with a broad peak coming down from 
the higher-temperature side. We can then see that the 
merged peak rapidly develops into a sharp peak associ- 
ated with the critical divergence at T/J ~ 0.8. 
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C^, and Ccross for the 3D Ising model of J' /J 



We next resolve the peak structure of the total specific 
heat using C,,, C±, and Ccmss- Figure 7(b) illustrates 
C||, which exhibits a broad peak for a small system size. 
Since J 3> J', this broad peak can be attributed to the 
spin fluctuation in the chain, which is governed by the 
energy scale J. As L increases, the peak temperature 
decreases to T 0.8 and the peak height itself increases 
in accordance with the criticality. On the other hand, 
C± in Fig. 7(c) shows a clear finite-size peak originating 
from the interchain fluctuation of the energy scale LJ'. 
As L increases, the peak temperature gradually increases 
from T/J ~ 0.4 to 0.8. Then, an interesting point about 
C± is that the shape of the peak is almost unchanged 
during shifting, in contrast to the 2D case where the 
peak considerably reduces its shape. Here, let us recall 
that, at sufficiently low temperatures, the effective spins 
are frozen in the chain direction form the 2D network. 
Thus, an important difference between the 2D and 3D 
cases is that, for 3D, the effective 2D Ising model in the 
small J' limit can involve the quasi-critical divergence, 
while for 2D, the specific heat of the effective ID Ising 
model shows no such divergence since there is no phase 
transition in the ID Ising model. In Fig. 7(d), we finally 
show Ccross, the peak of which develops near Tc and is 
smoothly connected to the critical divergence. 

In Fig. 8, we summarize the above size dependences of 
the peak temperatures for the specific heats. In the fig- 



ure, the horizontal axis indicates the scaled system size 
L~^/'^, where we have used v = 0.6301.^" We can see that 
the round peak of C„ comes from the higher temperature 
side, but it still does not reach the scaling region. On the 
other hand, we can see that Cj_ and Ccross are well fit- 
ted by linear functions, which are respectively shown as 
solid and broken lines in Fig. 8. The straightforward ex- 
trapolation yields Tc ~ 0.85, which is consistent with the 
precise estimation Tc — 0.834 based on the simulation up 
to the size 10 x 10 x 80 (the result is not presented here). 
A similar analysis of susceptibility was also reported in 
Ref. 8, where the peak temperature of the intrachain spin 
fluctuation behaves similarly to C±. The present result 
is consistent with the previous analysis of susceptibility.^ 
As can be seen Fig. 7(a), the divergences of Ccross and 
C± massively contribute to the critical divergence of the 
total specific heat C within a small system size. This 
suggests that Ccross Can be expected to be suitable for 
the finite-size analysis of the critical behavior as well, al- 
though C exhibits a rather complicated size dependence 
of the peak structure in the 3D case. 
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Fig. 8. Size dependences of the peak temperatures of C, Ch, C^, 
and Ccross for the 3D Ising model of J'/ J = 0.025. The scale of 
the horizontal axis follows L^^/" with v = 0.6301.'^'' 



4. Summary and Discussion 

We have studied the feature of the QID Ising model 
using Wang-Landau simulation. In order to treat the 
difference between the energy scales for the intra- and 
inter-chain directions, we have particularly introduced 
the two-dimensional energy space (e,,, ej^) corresponding 
to the intra- and inter-chain directions. We further de- 
composed the total specific heat C into contributions 
from intrachain fluctuation, d, interchain fluctuation, 
C_L , and the cross term of the intra- and inter-chain fluc- 
tuations, Ccross- Then the finite-size effect peculiar to 
the QID system is discussed on the basis of the two- 
dimensional DOS, and it was demonstrated that the 
chain flip conflguration plays an essential role for the low- 
temperature peak of specific heat. We have also analyzed 
the shift exponent of the peak of the specific heat, and 
then found that Ccross can capture the critical behavior 
more effectively than the total specific heat C, within 
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a relatively small system size. We have also discussed 
the qualitative difference between 2D and 3D cases in 
the low-temperature and small- J' limits; in the 2D case, 
the crossover of C± from an effective ID chain to a 2D 
model occurs rapidly at L ~ J /J'. In the 3D case, the 
effective 2D Ising model itself involves the quasi-critical 
divergence, because C± smoothly crossovers from the ef- 
fective 2D model into the 3D critical behavior. 

In this paper, we have analyzed the QID Ising model in 
the context of the two dimensional DOS. The actual com- 
putational cost for obtaining the two-dimensional DOS 
increases rapidly, with increasing system size, and then 
the cluster algorithm seems to be more efficient for a 
simulation of a larger system. However, the present de- 
scription based on the two-dimensional DOS provides the 
essential insight for the qualitative understanding of the 
low-energy excitations in the QID system. In addition, to 
suppress the finite-size effect peculiar to the QID system, 
the aspect ratio usually follows the ratio of anisotropic 
correlation length in analyzing the critical behavior. 
We can also see that a possible aspect ratio of the sys- 
tem is L' /L = {D — 1) J'/J, where D is the dimension of 
the system, L is the length of a chain, and L' is the size 
of the interchain direction. This is because the scale of 
the single-spin flip excitation and the chain flipped state 
can be on the same order 4J + 4{D-1)J' ~ A{D-1)J'L 
for J > J'. 
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